
#include "Lattice.h"
#include <string>
#include <sstream>

using namespace std;

#define Nx				64
#define Ny				64

#define Nu				0.002						//viscosidad (unidades LBM)
#define Tau				((6.0*Nu+1)/2.0)			//(S^{-1}_{nu} en taylor-green)
#define Omega			(1./Tau)					//factor de viscosidad de LBM (S_{nu} en taylor-green)
#define Rho				1.0							//densidad cte

#define U0				0.05
#define Kx				(2.0*M_PI/Nx)
#define Ky				(2.0*M_PI/Ny)
#define P0				0.0

#define Ntimes			149							//pasos temporales.
#define Nglobals		1							//periodo de escritura de globals.
#define Ngraph			(Nglobals)					//periodo de escritura de graficos (i.e. vorticidad).

void print_globals(Lattice lat) {
	cout << "RMS real velocity = " << lat.RMSvelocity() << endl;
	cout << "mean divergence = " << lat.divergence() << endl;
	cout << "energy = " << lat.energy() << endl;
	cout << "energy_k = " << lat.energy(FOURIER_SPACE) << endl;
	cout << "enstrophy = " << lat.enstrophy() << endl;
	cout << "palinstrophy = " << lat.palinstrophy() << endl;
}

void writeParameters(ostream& outStream) {
	outStream << "&INPUTS NX = " << Nx << ", NY = " << Ny << ", OMEGA = " << Omega << " /" << endl << endl;
	outStream << "Nu = " << Nu << endl;
	outStream << "Tau = " << Tau << endl;
	outStream << "Ntimes = " << Ntimes << endl;
	outStream << "Ngraph = " << Ngraph << endl;
}

int main(int argc, char *argv[]) {
	Lattice lattice;
	lattice.init(Nx, Ny, true, true, Rho, 1.0, 1.4, 1.7, 1.976284, 0); //se, sepsilon, sq, snu
	cout.precision(10);

	for(int i=0; i < Nx; i++)
		for(int j=0; j < Ny; j++) {
			double u0x = -U0*cos(Kx*i)*sin(Ky*j);
			double u0y = (Kx/Ky)*U0*cos(Ky*j)*sin(Kx*i);
			double P = (-0.25)*U0*U0*(cos(2.0*Kx*i) + pow(Kx/Ky,2.0)*cos(2.0*Ky*j)) + P0;
			double _rho = 1.0 + 3.0*P;
			vec<double> _u(u0x, u0y);
			double dxux = U0*Kx*sin(Kx*i)*sin(Ky*j);
			double dyux = -U0*Ky*cos(Kx*i)*cos(Ky*j);
			double dyuy = -U0*Kx*sin(Ky*j)*sin(Kx*i);
			double dxuy = U0*(Kx*Kx/Ky)*cos(Ky*j)*cos(Kx*i);
			lattice.setU(i, j, _u, dxux, dyux, dyuy, dxuy, _rho);
			//cout << _u.x << "\t\t" << lattice.getU(i,j).x << "\t\t" << _u.y << "\t\t" << lattice.getU(i,j).y << endl;
		}

// 	//print parameters
// 	cout << "--- Parameters ---" << endl;
// 	writeParameters(cout);
// 	cout << "--- Globals ---" << endl;
// 	print_globals(lattice);

	//LBM integration
// 	lattice.relaxRho();
// 	cout << endl << "--- Starting LBM ---" << endl;
	for(int n=0; n<Ntimes; n++) {
		if (n % int(Nglobals) == 0)
			cout << n << "\t" << lattice.energy()*4.0/(Nx*Ny*U0*U0) << endl;
		lattice.update();
	}
	return 0;
}
